Surrogate-Assisted Differential Evolution for the Design of Multimode Resonator Topology

The microstrip devices based on multimode resonators represent a class of electromagnetic microwave devices, promising use in tropospheric communication, radar, and navigation systems. The design of wideband bandpass filters, diplexers, and multiplexers with required frequency-selective properties, i.e., bandpass filters, is a complex problem, as electrodynamic modeling is a time-consuming and computationally intensive process. Various planar microstrip resonator topologies can be developed, differing in their topology type, and the search for high-quality structures with unique frequency-selective properties is an important research direction. In this study, we propose an approach for performing an automated search for multimode resonators’ conductor topology parameters using a combination of evolutionary computation approach and surrogate modeling. In particular, a variant of differential evolution optimizer is applied, and the model of the target function landscape is built using Gaussian processes. At every iteration of the algorithm, the model is used to search for new high-quality solutions. In addition, a general approach for target function formulation is presented and applied in the proposed approach. The experiments with two microwave filters have demonstrated that the proposed algorithm is capable of solving the problem of tuning two types of topologies, namely three-mode resonators and six-mode resonators, to the required parameters, and the application of surrogated-assisted algorithm has significantly improved overall performance.


Introduction
Microwave devices are widely used in many areas of communication and sensor networks.In particular, high-frequency selective devices based on microstrip resonators are characterized by a low cost in mass production and high-quality features [1].They can also be applied as sensors for measuring the electrical properties of materials [2,3].However, their main application is in the field of telecommunication, both tropospheric and space, as well as navigation [4].Different types of frequency-selective devices include bandpass filters, diplexers, multiplexers, and power dividers [5,6].
The traditional application field, which uses radio equipment, requires new promising designs, as the development of high-quality electronics relies on the continuous improvement of electrical characteristics and hence opens possibilities in size reduction [4,7].Therefore, studying the space of possible topologies of microstrip resonators is an important research field.The design features of the strip conductor allow for bringing together several lower oscillation modes in frequencies, which form a bandwidth [8].Based on the resonators it is possible to implement filters, including both bandpass and notch (stop band) filters [9,10].However, the microwave sensor design is also of particular interest [11,12].
The idea here is to use the changes in frequency response based on dielectric permittivity [13][14][15][16].These designs can be used to perform dielectric measurements of solids and liquids [17,18], as well as gas sensing [12] and noninvasive medical measurements [19,20].
The two-mode [1] and three-mode designs [21] are relatively simple, which, however, allow to perform frequency-selective tasks with high efficiency if properly tuned.Due to the small number of links in such devices, their size is relatively small, as well as the passband losses.The strip conductors usually have a rectangular form or a combination of rectangles.In a bandpass filter, the amount of coupling required between the resonances is created by etching one of the conductor corners in the shape of a square, or by adding a segment of a slotted line, or a strip coupling element, which are usually located at a 45 degree angle to the resonator axes [1].
The possible requirements for resonators could be different, for example, in some cases, it is important to cover a wide range of frequencies, and in other cases, it is important to obtain low losses.Due to their small size, microstrip resonators could be applied in small wireless networking devices, including radar sensors [22,23] and antenna sensors [2].However, it requires very accurate tuning and manufacturing of the devices-small changes in the size of the construction may significantly change the properties.Hence, here, we consider the problem of tuning the resonators to the required properties using software electrodynamic analysis tools.The capabilities of modern computing systems enable us to perform electrodynamic analysis many times while changing the geometric parameters of the microwave device.
The tuning and design of resonator topologies is a computationally expensive optimization problem due to the need to run simulations for every solution.There are various optimization techniques developed for expensive problems, for example, the surrogate approach, in which a model of the target function landscape is built [24].Based on the usage of Gaussian processes, which are also called kriging [25], efficient global optimization (EGO) was proposed [26].The EGO algorithm has found many applications in real-world, expensive optimization tools.
As the properties of the optimization problem to be solved are not known beforehand, and the derivative information is not available, searching for an optimal topology can be considered a black-box optimization task.A large variety of efficient black-box optimizers have been proposed in the area of evolutionary computation (EC) [27,28], including genetic algorithms (GAs) [29], particle swarm optimization (PSO) [30], and evolutionary strategies (ESs) [31].However, in the recent decade, differential evolution (DE) has become one of the most widely used numerical optimizers due to its simplicity and high efficiency [32].EC methods have also found many real-world applications [33].
The focus of this study is on combining the advantages of both evolutionary algorithms and surrogate-based optimization methods in order to solve a complex problem of multimode resonator design.In particular, a recently proposed variant of the differential evolution, L-SRDE, is combined with kriging in a local area around the best-known solution.The experiments are performed with two types of constructions, and the comparison results show that the proposed hybrid surrogate-assisted approach significantly outperforms both the differential evolution and the surrogated-based optimization.
The rest of this paper is organized as follows.The second section describes the microstrip resonators, their topologies, and their properties, as well as the DE and surrogatebased optimization algorithms.The third section contains a description of the proposed approach to automated microwave device tuning.The experimental set-up and results are shown in the third section, and the last section concludes this study.

Microstrip Resonators
Microstrip resonators are components of selection devices that are used in different types of communication devices and are mainly characterized by their AFC.In this study, planar resonators are considered, in particular, two constructions.The first one is U-shaped (slot-split rectangular), and the second one has U-shape with an additional cut (twice slot-split rectangular).The approach proposed here is not limited to these particular types of constructions and can be used for any other shapes.
Figure 1 shows the U-shaped microstrip resonator, in which the bandwidth is created by its two or three lowest oscillation modes [34,35].The width and depth of the conductor cut influence the position of the lower natural frequencies, and tuning them allows bringing them closer.The bandwidth of a resonator can be determined as follows: where ΔF = F h − F l , and F h and F l are the high-frequency and low-frequency boundaries of the bandwidth, measured at a −3 dB level from minimal losses, and F 0 is the center frequency of the bandwidth.Some of the previous studies on three-mode resonators of such construction have shown that the minimum bandwidth is around 45%, and the maximum is 71% when F 0 = 3 GHz [35].Such structure is widely used due to its simplicity, and it finds application in various types of devices, such as sensors [3,36,37].The results of the numerical electrodynamic modeling of such a resonator on a TBNS dielectric plate (thickness h = 1 mm, relative permittivity ε = 80) are very close to the parameters of the manufactured prototype [35].Figure 1 also shows the geometric characteristics of the U-shaped 3-mode resonator.In particular, there are 13 parameters, set to the following values: 1.
x 3 = 2.20 mm-width of the center rectangle;
x 5 = 3.50 mm-length of right side port connection (always equal to x 1 ); 6.
y 1 = 0.98-relative height of the left side port connection; 7.
y 3 = 2.80 mm-distance from the bottom edge of the dielectric plate; 9.
y 4 = 18.90 mm-height of the left rectangle; 10. y 5 = 12.80 mm-height of the center rectangle; 11. y 6 = 18.90 mm-height of the right rectangle; 12. y 7 = 0.98-relative height of the right side port connection; 13. y 8 = 2.80 mm-distance to the top edge of the dielectric plate.
It is worth mentioning that unlike previous work [38], where the left-and right-hand sides of the resonator were symmetric and mirroring was applied, here, we allow asymmetric geometry.This leads to more parameters to be set but allows wider investigation of the proposed approach, as it can be further applied to design diplexers and multiplexers.From the 13 parameters listed above, only 8 are tunable, which are the rectangles parameters x 2 , x 3 , x 4 , y 4 , y 5 , y 6 and port positions y 1 and y 7 .The parameters described here result in a resonator with F 0 = 1.498GHz and ΔF = 0.618, as shown in Figure 2. The second resonator considered in this study is shown in Figure 3, and its AFC is presented in Figure 4.This is a twice slot-split rectangular microstrip resonator, capable of 6-mode operation.In this case, the theoretical investigations are in agreement with the characteristics, measured on a manufactured prototype, meaning that the experiments could be performed automatically using only electrodynamic modeling tools.Previous experiments with this type of structure have shown that at F 0 = 3 GHz, it allows one to reach a bandwidth of around 80%.
All 17 parameters of the second structure, their values, and meanings are listed below: 1.
y 1 = 0.98-relative height of the left-side port connection; 9.
y 2 = 0.02 mm-width of the port connections; 10. y 3 = 2.40 mm-distance from the bottom edge of the dielectric plate; 11. y 4 = 28.00mm-height of the first rectangle; 12. y 5 = 21.40 mm-height of the second rectangle; 13. y 6 = 15.70 mm-height of the third rectangle; 14. y 7 = 21.40 mm-height of the fourth rectangle; 15. y 8 = 28.00mm-height of the fifth rectangle; 16. y 9 = 0.98-relative height of the right side port connection; 17. y 10 = 3.00 mm-distance to the top edge of the dielectric plate.
From the 17 parameters listed above, only 12 are tunable, which are the rectangular parameters x 2 , x 3 , x 4 , x 5 , x 6 , y 4 , y 5 , y 6 , y 7 , y 8 and port positions y 1 and y 7 .The parameters described here result in a resonator with F 0 = 1.832GHz and ΔF = 0.957, as shown in Figure 4.The two types of resonator structures are used in this study for experiments with automatic topology tuning.The purpose of the search process will be to reach a certain frequency range and bandwidth.It was mentioned that for the second example, the twice slot-split rectangular microstrip sensor, it is possible to reach 6-mode operation, but the initial solution will have only 4 modes, as shown in Figure 4. So, one of the goals of the optimization method is to find the desired number of modes.The solution evaluation method is described in further sections.

Continuous Optimization Methods
The topology of the constructions described above is controlled by a number of numerical values, and tuning them allows for achieving different frequency-selective properties.However, the evaluation of a particular set of parameter values is expensive, i.e., time-consuming, so an efficient optimization technique should be applied in order to obtain the constructions with desired properties.
In general, an optimization problem can be formulated as follows: given a target function f (x), find an optimal vector x opt so that x ∈ R D , where D is the search space dimension, and x j ∈ [x j,min , x j,max ] with x j,min and x j,max are the lower and upper boundaries.
The field of optimization methods nowadays is vast and includes a variety of approaches, depending on the problem.In the case considered in this study, the derivative information is not available, so a zero-order method that requires only target function values should be used.The existing derivative-free optimization techniques could be divided into local and global search methods; the former may be more cost-efficient but often becomes stuck in local optimum, whereas the latter have better exploration capabilities at the cost of additional target function evaluations.An efficient optimization tool should combine the advantages of both of these approaches.

Surrogate-Based Optimization
Solving a real-world engineering problem usually requires building a mathematical model of the considered system.If the theoretical knowledge about the system is available and usable, then the model can be built using it; however, in many cases, this is not applicable.When dealing with complex optimization problems, the only possible way of building a mathematical model could be through experimental data.The earliest studies have relied on the design of the experimental approach [39,40], but a more recent approach is to apply Bayesian methods, such as Gaussian processes (GPs), also known as kriging [41].Based on the idea of applying GPs, the Efficient Global Optimization (EGO) algorithm was proposed [26].In EGO, the built model is used as a surrogate, i.e., the optimum is found using it, and it becomes the new point to be evaluated.As kriging allows evaluating not only the mean µ but also the standard deviation σ in each point, in EGO, it is possible to use several types of criteria to determine the next point: 1.
Expected Improvement (EI), using the function i.e., combining the information about the best point, µ and σ.
Here, Y is a random variable with normal distribution and parameters µ and σ, f min is the best-known function value, and E is the expectation.One of the disadvantages of EGO is that generating each new point requires solving another optimization problem using the surrogate model, and training the surrogate on the available data becomes more and more time-consuming as the dataset size and problem dimension grow.Nevertheless, the EGO approach has found many applications and proved to be an efficient search method for complex problems [42].

Differential Evolution
The evolutionary algorithms are optimization techniques, inspired by the process of natural selection.The earliest versions of EAs, such as genetic algorithms, have been shown to perform well on a variety of global optimization problems in different domains, including binary, combinatorial, integer, and numeric optimization [43].There have been many other classes of methods proposed besides GAs, for example, evolutionary strategies (ES) and differential evolution (DE) [44].The latter today is one of the most efficient and widely used global numeric optimizers due to simplicity and high efficiency [45].In the latest optimization competitions of the Congress on Evolutionary Computation (CEC), the DE-based approaches occupy not just the leading positions but most of the submitted algorithms [46,47].
Differential evolution is a population-based optimizer, which means that it uses a set of N randomly generated vectors x i,j , i = 1, ..., N, j = 1, ..., D, where D is the problem dimension.The main operation in DE is the difference-based mutation, one of the most widely used strategies is called current-to-pbest/1 [48]: (3) where x i is the i-th vector in the population of N individuals, and each vector consists of D variables; v i is called the donor vector; r1 and r2 are random indices from the population; pbest is one of the p% best individuals; and F is called the scaling factor.This mutation consists of two parts, exploitation, i.e., moving toward one of the p% best solutions with the first difference, and exploitation, that is, moving in a random direction with the second difference.This allows one to use both trends in DE.
Mutation is followed by crossover; the most popular method is the binomial crossover, implemented as follows: where u i is the trial vector, jrand is required to make sure that at least one component is taken from the donor vector, and Cr is the crossover rate parameter.The generated trial vectors are also checked to be within the bounds and corrected if needed [49].
The trial solutions u i are then evaluated and compared with the target ones in the selection step: That is, if the newly generated solution is better, it replaces the old one.
All modern DE variants rely on one of the parameter adaptation techniques.The reason for this is that the three main parameters, population size N, scaling factor F, and crossover rate Cr, significantly influence performance [50].As the DE is designed to be a general-purpose black-box optimizer, the best way to set the parameters is to determine them based on how well certain values perform during the search.For the last 10 years, one of the most efficient parameter adaptation techniques for DE is the successful history-based adaptation (SHA), proposed in [51].The L-SHADE algorithm [52], which tunes F and Cr with SHA and linearly decreases the population size, is a baseline approach for many modern DE variants [46].
Recently, there have been several attempts to change the general scheme of L-SHADE, for example, in the L-NTADE algorithm [47], where two populations are used.The successful history-based adaptation was shown to be biased [53], and a hyperheuristic approach was applied to design the parameter adaptation method automatically [54].The algorithm proposed in this study is based on these results.

Previous Work on Microstrip Multimode Resonator Design
In our previous study [38], the same problem of microstrip multimode resonator topology design was considered.In that work, three types of constructions were considered: 3-mode, 5-mode, and 6-mode resonators.The main task was the same: given a baseline topology, tune its bandpass to the desired frequency range while keeping the same width.
To optimize the topology, several approaches were considered, namely simulated annealing (SA), continuous genetic algorithm (CGA), and a hybrid between the genetic algorithm and differential evolution with parameter adaptation.A special fitness function was designed to evaluate the solutions, combining information about the number of modes, reflection losses, and the operating frequency range.Also, a heuristic approach was adopted to obtain an approximation of the desired solution given the initial and desired frequency-this was performed by scaling the whole topology.
The drawbacks of that study were that, in fact, the local search was performed around the scaled solution, and the mesh was fixed during the simulations.The computational resource was limited to 500 function evaluations.That is, although it was shown that it is possible to tune resonators to desired characteristics with this approach, it is not general.For example, if instead of tuning to a different frequency the aim is to increase the width of the passband, the scaling would not help, and this would require solving a global optimization problem.The methods used to overcome these limitations are described in the next section.

Proposed Approach: Target Function Formulation
The main focus of this story is on the development of an optimization technique, which would be general and capable of tuning different types of microwave structures.To achieve this, first of all, we determine the target function, also known as the fitness function in the evolutionary computation community.The target function will be used to determine the quality of the solution, i.e., its difference from the desired one, with a single numeric value.In the previous study, the quality of the solutions was determined with five different metrics, which were then combined into a single one [38].Here, we follow a similar approach, but the metrics are changed to have the same scale-our preliminary experiments showed that this is crucial for obtaining high-quality results.
Let us suppose that the goal is to tune the resonator so that it would have the following two values determining the bandwidth: 1.
High-frequency F des u .where "des" stands for "desired".Hence, the width can be calculated as follows: The denominator is also called the center frequency, F des 0 = 0.5(F des u + F des l ), and the nominator is ΔF des = F des u − F des l .The same set of values is calculated for the actual topology based on the results of electrodynamic modeling.These actual values are F act l , F act u , F act 0 , ΔF act , dF act .The optimization has the following five objectives: 1.
The desired and actual center frequencies of bandpass should be the same, The desired and actual bandwidth should be the same, dF des = dF act ; 3.
The number of actual modes should be greater or equal to the desired ones, N M act ≥ N M des ; 4.
The frequency dependence of losses in reflection should not exceed the threshold level within the bandwidth, i.e., the actual losses threshold should be LT act < LT des ; 5.
The width at the threshold level LT des should be greater than the scaled desired width, i.e., Here, F act u,T and F act l,T are the frequencies where the reflection losses hit the LT des level, and CT ∈ (0, 1] determines the allowed narrowness at this level compared with dF des .The LT des was set to −14 dB in all experiments.For the last objective, we determine Note that the three last objectives could also be formulated as constraints.The mathematical formulation of all five objectives is given below.

1.
f it < 1, otherwise 0; Also note that in the fourth objective, both LT act and LT des are negative, and, for example, LT act = −13 is penalized.Finally, the fitness of a solution is just a sum of all five: An example of calculating all of these values will be given in the experimental section.

Proposed Approach: Differential Evolution Modification
As mentioned above, in this study, a combination of the evolutionary algorithm, namely a variant of differential evolution, and surrogate-based optimization is proposed.First, the used DE algorithm is described, and then its hybrid with kriging is as well.
The DE algorithm proposed here is based on several works on parameter adaptation techniques for differential evolution.In particular, as DE is highly sensitive to setting the scaling factor F, in [53], it was shown that the existing successful history-based adaptation can be modified by biasing generated F to higher values.However, it is possible to design adaptation techniques automatically using genetic programming, as shown in [54].Based on the findings of the paper in [55,56], a hyperheuristic approach to automatically design parameter adaptation techniques was implemented with a simple rule for tuning scaling factor F based on success rate.
The success rate SR is determined by the ratio of the number of improved solutions NS to the current population size N: The number of improved solutions NS is calculated as the number of times when the selection (Equation ( 5)) step was successful.In particular, the mean for sampling scaling factor F values can be calculated as the c-th order root of the success rate [57]: where c = 4 is a typical setting of the parameter.Finally, the sampling for F is performed as follows: where randc(m, s) is a Cauchy-distributed random number with location parameter m and scale parameter s.If the sampled F value is smaller than 0, it is sampled again, and if it is larger than 1, it is set to 1.The proposed algorithm called SRDE (Success-Rate-based Differential Evolution) also uses the external archive, the same as in the SHADE [51] algorithm.The first step is the initialization, which is performed based on the given baseline solution b j (parameters from Section 2.1 are used here).In particular, the population x i,j and the archive a i,j (i = 1, 2, ...N, j = 1, 2, ...D) are generated with normal distribution as follows: x i,j = randn(b j , 0.05 • (xmax j − xmin j )) (13) where randn(m, s) is a normally distributed random variable with mean m and standard deviation s, and xmax j , xmin j are the upper and lower bounds for variable j.In other words, the scaled range 0.05 • (xmax j − xmin j ) is used as the standard deviation.The archive a i,j is sampled in the same way.The mutation strategy relied on selective pressure, which changes the probabilities of individuals being chosen for mutation based on their fitness, as described in [58].In particular, the exponential selective pressure was applied to the r1 index in Equation (3).The rank of a particular individual in an array sorted by fitness is calculated as follows: where kp is the selective pressure parameter.The resulting probabilities are calculated as After initialization, the main cycle begins.F values are generated as described above, and the crossover rate parameter is fixed as Cr = 0.9.The mutation strategy used is the current-to-pbest with archive, the same as in Equation ( 3), but the last vector x r2 can be chosen from the joined set of population and archive.After binomial crossover (Equation ( 4)), the solutions are checked to be within the search bounds using the midpointtarget method: During the selection step, the solutions, which are replaced in the population, are put into the archive, replacing randomly chosen ones.Also, the number of successful replacements is calculated, as well as the success rate SR.This completes the single generation (iteration) of SRDE.

Proposed Approach: Surrogate-Assisted SRDE
One of the goals of the current study is to combine the fast convergence of the DE algorithm with the explorative properties of surrogate-based optimization.In particular, the hybrid approach, called surrogate-assisted SRDE (SA-SRDE), is proposed.In this algorithm, the basic algorithm is the SRDE, but in addition, the kriging model is built using the evaluated points, and every generation, one of the worst points in the SRDE population is replaced by the point found on the surrogate model.
The main steps of the proposed algorithm can be described as follows: 1.
Train a surrogate model KRG([x ds i , y ds i ]) on the dataset; 6.
Apply acquisition function using best known solution x b to obtain x a f = AF(KRG, x b ) and evaluate fitness y a f = f it(x a f ); 7.
Replace the worst individual x w in the population with x a f if y a f < y w ; 8.
If the termination condition is not met, go to step 2; 9.
Return best solution x b , y b .
The size of the dataset M may differ depending on the number of improved points.
Step 6 requires additional explanation, as it contains the acquisition function, which is another SRDE algorithm applied to the surrogate model.The search range is set around the best-known solution x input b . The steps of AF(KRG, x input b ) are described below: 1.
Use the best-known solution to set the temporary search boundaries: Initialize the SRDE algorithm (Equation ( 13)) and evaluate fitness using the LCB approach (Section 2.3): Generate mutant vectors v i (Equation 3)) and trial vectors u i (Equation ( 4)), apply bound constraints (Equation ( 15)); 4.
If the termination condition is not met, go to step 2; 6.
Return best solution x b , f it(x b ).
In Figure 5, the flowchart of the proposed approach is shown.If the surrogate is not used (SRDE), only the left part of the scheme is active; otherwise, both are active (SA-SRDE).Instead of the LCB approach, other variants, such as SBO or EI, could be applied.Note that the EGO algorithm uses the SLSQP approach to perform optimization on the surrogate model, but a combination of LCB and SRDE has shown better results in preliminary experiments.
In the next section, the details of the performed computational experiments are given, including parameter settings of the used algorithms, as well as the final results, approach comparison, and the applicability of the results for constructing various devices.

Results and Discussion
The number of experiments in this study is mainly limited by the computational complexity of the electrodynamic modeling of the resonators.Evaluating a single set of parameters may take several minutes, so finding a topology with acceptable characteristics may take several days.The optimization algorithm was developed in Python, and the Surrogate Modeling Toolbox (SMT) [59,60] was applied for EGO and kriging.The experiments ran on six computers with Intel Core i9 13900 processors and 32Gb RAM.

Experimental Set-Up
There were two topologies considered (see Section 2.1), and several experiments were performed for each of them.In particular, the following tests were conducted: 1.
EGO algorithm on the first topology with recommended parameters from SMT; 2.
Standard SRDE algorithm on the first topology; 3.
SA-SRDE algorithm with LCB-based acquisition function on the first topology; 4.
Standard SRDE algorithm on the second topology; 5.
SA-SRDE algorithm with LCB-based acquisition function on the second topology; For the first structure, the following required parameters were set: F des l = 1.5 GHz and F des u = 3.0 GHz, and for the second structure, F des l = 1.2 GHz and F des u = 3.2 GHz.During every experiment, all the evaluated parameters were stored for later analysis, as well as the time intervals to evaluate solutions and run the algorithm.The search ranges for every variable for the first and second topologies, as well as the baseline solution xbas parameters, are shown in Tables 1 and 2. 0.00 0.98 1.00

Results on the First Structure
Tables 1 and 2 only show the parameters, which were actually tuned during the optimization process.The ranges were set with enough margin so that the search algorithm would not get stuck because of them.
The EGO algorithm was used with the following settings: initial number of points: 25; total number of target function evaluations (simulations): NFE max = 1000; and criterion: expected improvement (EI).The parameters of the SRDE and SA-SRDE shared the same parameters regarding the differential evolution part, and they are listed below: 1.
Computational resource NFE max = 1000 evaluations for first structure, NFE max = 1500 for second structure; 2.
Number of independent runs NR = 15 for first structure, NR = 10 for second structure; 3.
Population size N = 25 for first structure, N = 36 for second structure; 4.
The best target function values for each of the 15 runs in the first three experiments are shown in Table 3.As can be seen from Table 3, the EGO algorithm was not able to achieve good results in most of the runs; however, both SRDE and SA-SRDE performed much better.SA-SRDE achieved the best mean, standard deviation, median, and maximum values but obtained the same result as SRDE in terms of the best found solution out of 15 runs.The best values in Table 3 are shown in bold.Table 4 contains the statistical analysis of the results, shown above.The Student's t-test and Mann-Whitney U test show that both SRDE and SA-SRDE outperform EGO; however, the comparison between SRDE and SA-SRDE does not show any statistically significant differences between these two.Nevertheless, as Table 3 has shown, SA-SRDE is better on average.
Figure 6 shows the convergence graphs of the three tested algorithms.The wide lines are the average over all runs.As Figure 6 demonstrates, the EGO algorithm performs much worse than the SRDE and SA-SRDE.Comparing SRDE and SA-SRDE, it can be noted that after 300 evaluations, the SA-SRDE starts to show better results on average.

Results on the Second Structure
The results of the second set of experiments with the second structure are shown in Tables 5 and 6 contains the statistical comparison.The experiments with EGO on the second structure were not performed, because this method showed poor results in the first set of tests.The best values in Tables 5 and 6 are shown in bold.As Table 5 shows, the SA-SRDE outperformed the baseline approach in most of the basic statistical characteristics.The statistical tests also show very small p-values, around 8.7 × 10 −2 , meaning that the difference is rather significant.
Figure 7 is similar to Figure 6 and shows the convergence graphs of the two tested algorithms for the second structure.The wide lines are the average over all runs.Figure 7 demonstrates that SA-SRDE outperforms SRDE after 500 evaluations and achieves better results on average at the end.

Time Requirements Analysis
Considering the analysis above, we may conclude that the SA-SRDE performs better than the baseline approach SRDE due to the usage of the surrogate approach.However, it is important to consider that the SA-SRDE requires additional computational effort in order to build models and optimize them.The time requirements in seconds of these two methods are shown in Tables 7 and 8.The comparison of time, required to finish the optimization in Tables 7 and 8, demonstrates that the SA-SRDE does not always work longer.In particular, for the first structure, the SA-SRDE worked around 10% longer than the baseline algorithm SRDE on average, but for a more complex second structure, the times of SA-SRDE are smaller.This can be explained by the fact that the tuning is performed towards higher frequency ranges, and it results in smaller structures in terms of geometric characteristics.Smaller structures tend to require less computational resources to be modeled, so as SA-SRDE performs better and replaces worse solutions, it generally evaluates smaller structures more often, resulting in an advantage in computation time, as Table 8 shows.However, this advantage would disappear if the goal would be to tune to lower frequency ranges.

Analysis of Best Found Solutions
To show the quality of the best found solutions in this section the best solutions found for both structures are shown.Figure 6 shows the first structure, optimized with SA-SRDE, 15-th run, and Figure 6 shows its AFC.Table 9 contains the geometric parameters.Figure 5 shows that the optimization process results in an asymmetric structure of the resonator.The proposed approach allows one to obtain full symmetry, but it requires increased computational resources.The AFC characteristics of the found solution are almost ideal: the frequencies are very close to the desired ones, and the S 11 curve does not reach the −14 dB level.
Figures 8 and 9 show the geometry and AFC of the optimized second structure from run 2 of SA-SRDE, and Table 10 contains the geometric parameters.Unlike the first structure, here, in Figure 10, it can be seen that the topology is almost perfectly symmetric, and the values in Table 10 support this statement.As for the AFC in Figure 11, unlike the one shown in Figure 4, here, we may clearly observe six modes, and the frequency range is even larger than the required 1.2 and 3.2 GHz.It shows that there is a potential room for further improvement of this variant.

Conclusions
In this paper, a complex engineering optimization problem of tuning the conductor's topologies of the microstrip multimode resonators was considered.This task is characterized by high computational complexity, as it requires electrodynamic modeling of the microwave devices, which takes minutes of time even on modern hardware.To solve this problem a hybrid surrogate-assisted differential evolution algorithm was proposed and tested on the problem of tuning two resonators with different topologies.The experiments have shown that this approach outperforms the baseline algorithm without the kriging model and allows one to find the efficient topologies of microwave devices.
It is important to mention that the approach proposed in this study is not limited to the design of microstrip resonators and can be utilized in other areas of complex engineering design, where the numeric parameters should be tuned.In the area of microwave devices, the algorithm can be applied to design diplexers, multiplexers, or power dividers, and the only things that should be changed are the search ranges and fitness calculation.As for developing more efficient search algorithms, unfortunately, it is difficult to search for more efficient optimizers for this type of expensive problem, so further developments should be performed on simpler problems that have similar landscape properties to the problems considered in this study.Determining these properties could be one of the directions of further studies.

21 Figure 11 .
Figure 11.AFC of the optimized U-shaped microstrip resonator with additional cut.

Table 1 .
Search range and baseline solution for the first structure, in mm.

Table 2 .
Search range and baseline solution for the second structure, in mm.

Table 3 .
Results of experiments with the first structure.

Table 4 .
Statistical comparison of first three experiments.

Table 5 .
Results of experiment with the second structure.

Table 6 .
Statistical comparison of the two last experiments.

Table 7 .
Time requirements of first three experiments in seconds.

Table 8 .
Time requirements of the last two experiments in seconds.

Table 9 .
Optimized parameters of the first structure in mm.

Table 10 .
Optimized parameters of the second structure in mm.